home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
EnigmA Amiga Run 1999 March
/
EnigmA AMIGA RUN 35 (1999)(G.R. Edizioni)(IT)[!][issue 1999-03].iso
/
earcd
/
devel
/
vbcc-68k-src
/
machines
/
amiga68k
/
libsrc
/
math
/
math_040
/
sinh.s
< prev
next >
Wrap
Text File
|
1999-01-01
|
5KB
|
120 lines
*
* $VER: sinh.s 33.1 (22.1.97)
*
* hyperbolic sine
*
* Version history:
*
* 33.1 22.1.97 (c) Motorola
*
* - snipped from M68060SP sources
*
machine 68040
fpu 1
XDEF _sinh
XDEF @sinh
XREF @exp
XREF @expm1
*************************************************************************
* sin(): computes the hyperbolic sine of a normalized input *
* *
* INPUT *************************************************************** *
* fp0 = extended precision input *
* *
* OUTPUT ************************************************************** *
* fp0 = sinh(X) *
* *
* ACCURACY and MONOTONICITY ******************************************* *
* The returned result is within 3 ulps in 64 significant bit, *
* i.e. within 0.5001 ulp to 53 bits if the result is subsequently *
* rounded to double precision. The result is provably monotonic *
* in double precision. *
* *
* ALGORITHM *********************************************************** *
* *
* SINH *
* 1. If |X| > 16380 log2, go to 3. *
* *
* 2. (|X| <= 16380 log2) Sinh(X) is obtained by the formula *
* y = |X|, sgn = sign(X), and z = expm1(Y), *
* sinh(X) = sgn*(1/2)*( z + z/(1+z) ). *
* Exit. *
* *
* 3. If |X| > 16480 log2, go to 5. *
* *
* 4. (16380 log2 < |X| <= 16480 log2) *
* sinh(X) = sign(X) * exp(|X|)/2. *
* However, invoking exp(|X|) may cause premature overflow. *
* Thus, we calculate sinh(X) as follows: *
* Y := |X| *
* sgn := sign(X) *
* sgnFact := sgn * 2**(16380) *
* Y' := Y - 16381 log2 *
* sinh(X) := sgnFact * exp(Y'). *
* Exit. *
* *
* 5. (|X| > 16480 log2) sinh(X) must overflow. Return *
* sign(X)*Huge*Huge to generate overflow and an infinity with *
* the appropriate sign. Huge is the largest finite number in *
* extended format. Exit. *
* *
*************************************************************************
T1 dc.l $40C62D38,$D3D64634 ; 16381 LOG2 LEAD
T2 dc.l $3D6F90AE,$B1E75CC7 ; 16381 LOG2 TRAIL
_sinh
fmove.d (4,sp),fp0
@sinh
fmove.x fp0,-(sp)
move.l (sp),d1
move.w (4,sp),d1
move.l d1,a1 ; save (compacted) operand
lea (12,sp),sp
and.l #$7FFFFFFF,d1
cmp.l #$400CB167,d1
bgt.b .SINHBIG
;--THIS IS THE USUAL CASE, |X| < 16380 LOG2
;--Y = |X|, Z = EXPM1(Y), SINH(X) = SIGN(X)*(1/2)*( Z + Z/(1+Z) )
fabs.x fp0
move.l a1,-(sp) ; {a1}
jsr @expm1 ; FP0 IS Z = EXPM1(Y)
move.l (sp)+,a1 ; {a1}
fmove.x fp0,fp1
fadd.s #$3F800000,fp1 ; 1+Z
fmove.x fp0,-(sp)
fdiv.x fp1,fp0 ; Z/(1+Z)
move.l a1,d1
and.l #$80000000,d1
or.l #$3F000000,d1
fadd.x (sp)+,fp0
move.l d1,-(sp)
fmul.s (sp)+,fp0 ; last fp inst - possible exceptions set
rts
.SINHBIG
cmp.l #$400CB2B3,d1
bgt .t_ovfl
fabs.x fp0
fsub.d (T1,pc),fp0 ; (|X|-16381LOG2_LEAD)
clr.l -(sp)
move.l #$80000000,-(sp)
move.l a1,d1
and.l #$80000000,d1
or.l #$7FFB0000,d1
move.l d1,-(sp) ; EXTENDED FMT
fsub.d (T2,pc),fp0 ; |X| - 16381 LOG2, ACCURATE
jsr @exp
fmul.x (sp)+,fp0 ; possible exception
.t_ovfl
rts